from pylab import *

d   = loadtxt('../outputs/xe_rhs.dat')
z   = d[:,0]
lna = d[:,1]
xe  = d[:,2]
rhs = d[:,3]

# Redshifts are outputed each time ODE_rec is called. 
# If the GSL ODE solver tries to use a smaller stepsize h, then the outputed redshifts {z_i}
# are not neccessarily in monotonically decreasing, and that's why I use numpy.argsort() to
# obtain idx_sorted from {z_i}
idx_sorted = argsort(z)

fig = figure(figsize=(12,6))

subplot(1,2,1)
plot(z[idx_sorted], rhs[idx_sorted], 'b-', label='$RHS$')
plot(z[idx_sorted], xe[idx_sorted], 'r-', label='$Xe$')

ax = gca()
ax.set_xlim(z.max(), z.min())
ax.set_xlabel(r'$z$',fontsize=14)
ax.set_ylabel(r'$\frac{d X_e}{d \ln a}$',fontsize=14)
ax.grid()

legend(loc='best', fontsize=14)


subplot(1,2,2)
plot(lna[idx_sorted], rhs[idx_sorted], 'b-', label='$RHS$')
plot(lna[idx_sorted], xe[idx_sorted], 'r-', label='$Xe$')

ax = gca()
ax.set_xlabel(r'$\ln a$',fontsize=14)
#ax.set_ylabel(r'$\frac{d X_e}{d \ln a}$',fontsize=14)
ax.grid()

legend(loc='best', fontsize=14)

tight_layout()

show()
